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g Abstract 

'^ First principles based phase diagram calculations were performed for the octahedral-interstitial 

O solid solution system aZrOjf {aZr[ ]i^xOx', [ ]=Vacancy; < X < 1/2). The cluster expansion 

method was used to do a ground state analysis, and to calculate the phase diagram. The predicted 

J> diagram has four ordered ground-states in the range < X < 1/2, but one of these, at X=5/12, 

o 

^^ is predicted to disproportionate at T~ 20K, well below the experimentally investigated range 

m 

■^ Tf« 420K. Thus, at T ^ 420K, the first-principles based calculation predicts three ordered phases 

l> 

(^ rather than the four that have been reported by experimentalists. 

> 

X 

S3 



Key words: ZtOx', 'ZiT suboxides; Zircalloy; First Principles; Phase diagram calcula- 
tion; vacancy-interstitial ordering; order-disorder; alloy theory. 



I. INTRODUCTION 

Zircalloy is used as nuclear fuel-rod cladding in light water reactors, but it is metastable 
with respect to oxidation by the UO2 fuel.'^'Sl 

Oxidation of zircalloy transforms it from the high-temperature (high-T), oxygen-poor, bcc 
solution (/3Zrx) into the low-T, oxygen-rich, hep-based solution (aZrOx)- At temperatures 
between about 1173K and 573K various ordered phases have been reported.'^'^ 

Octahedral interstitial ordering of oxygen (O), and vacancies ([ ]) in aZrOx 
{aZr[ ]i^xOx, < X < 1/2) increases microhardness^ and brittleness^, and therefore, 
promotes stress corrosion cracking. Order-disorder transitions were studied via heat capac- 
ity measurements: Aral and HirabayashP^ studied alloys with 0/Zr ratios of 0.16 and 0.24 
at 473K-973K; Tsuji and AmaysP^ made similar measurements on alloys with 0/Zr ratios 
of 0.0, 0.10, 0.13, and 0.24, at 325K-905K. 

Aral and HirabayashP^ achieved a high degree of long-range ordering in samples that 
were cooled from 623K to 523K, during a period of about one month, which indicates a high 
mobility of oxygen in aZr[ ]i^xOx, even at such modest temperatures; hence a system 
that is highly reactive even at such moderate temperatures. 

A recent computational study ^^ reported ground-state structures and order-disorder tran- 
sition temperatures for ZrgO and ZraO, but did not present a calculated phase diagram, or 
report if the calculated order- disorder transitions in ZrgO and Zt^O are first-order, as ex- 
periment indicates, or continuous. 

The results presented below are mostly consistent with experimental studies with respect 
to the comparison between computationally predicted ground-state (GS) structures and 
reported (experimental) low-temperature (T ^ 500-ft") ordered phases, with the exception 
that in the range 0.25 ^ X ^ 0.5 the calculations predict only two ordered phases at T150i^, 



rather than the three called a'g) c^s ^^^ ^'i i^ Aral and Hirabayashi (1976).'^ 

Experimental values for the maximum solubility of O in Zr, X^ax; range from: X^ax ~ 
29 at. V^ to X^ax ~ 35 at. 9^ and X^ax -40 at. 7c!^" ' ^^ ' " l. The first-principles results 
presented here support a higher value; i.e. X^ax > 1/2. This may reflect an insufficiently 
negative calculated value for the formation energy of monoclinic Zr02. 

II. METHODOLOGY 

A. Total Energy Calculations 



Formation energies, AiJj (Fig. uj) were calculated for fully relaxed hep aZr, hep aZrO 
(hep aZr with all octahedral interstices occupied by O), and 96 aZr[ ]i-nO„ supercells 
of intermediate composition. All calculations were performed with the density functional 
theory (DFT) based Vienna ah initio simulation program (VASP, version 449^^'^ using 
projector- augmented plane-wave pseudopotentials, and the generalized gradient approxima- 
tion for exchange and correlation energies. Electronic degrees of freedom were optimized 
with a conjugate gradient algorithm, and both cell constant and ionic positions were fully 
relaxed. Pseudopotential valence electron configurations were: Zr^^: 4s4p5s4d; O^: Ss^Sp"^. 

Total energy calculations were converged with respect to k-point meshes by increasing 
the density of k-points for each structure until convergence. A 500 eV energy cutoff was 
used, in the "high precision" option which guarantees that absolute energies are converged 
to within a few meV/site (a few tenths of a kJ/site of exchangeable species; O, [ ]). Residual 
forces were typically 0.02 eV or less. 

Calculated formation energies, A£'j, relative to a mechanical mixture of aZr + aZrO, 
for the 96 aZr[ ]i^nOn supercells are plotted as solid circles in Fig. Il| Values of AEj are, 

/S.Ef = {Estr — EaZr — EaZro) / C^) (l) 

where: Estr is the total energy of the aZr[ ]i_nO„ supercell; EaZr is the energy/atom of 
aZr] EaZrO is the energy/atom of aZrO. 
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FIG. 1: Comparison of VASP (solid circles) and CE (larger open squares, red online) formation 
energies, AE'j, and a ground-state analysis on structures with 18 or fewer octahedral-interstitial 
sites (smaller open squares, blue online). Extension of the convex hull towards the formation 
energy of monoclinic zirconia, Zr02, indicates that the four ordered GS at X=l/6, 1/3, 5/12 and 
1/2 are also predicted to be GS of the Zr-0 binary. 



B. The Cluster Expansion Hamiltonian 

The cluster expansion, CEP^, is a compact representation of the configurational total 
energy. In the aZr[ ]i-xOx system, the solid solution configuration is described by pseu- 
dospin occupation variables a,, which take values cij = — 1 when site-i is occupied by [ ] 
and cTj = +1 when site-i is occupied by O. 

The CE parameterizes the configurational energy, per exchangeable cation, as a polyno- 
mial in pseudospin occupation variables: 



Eia) = J2^iJAll^^) (2) 

i \iei' I 

Cluster £ is defined as a set of lattice sites. The sum is taken over all clusters i that are 
not symmetrically equivalent in the high-T structure space group, and the average is taken 
over all clusters C.' that are symmetrically equivalent to i. Coefficients Ji are called effective 
cluster interactions, ECI, and the multiplicity of a cluster, m^, is the number of symmetrically 
equivalent clusters, divided by the number of cation sites. The ECI are obtained by fitting 
a set of VASP FP calculated structure energies, {Estr}- The resulting CE can be improved 
as necessary by increasing the number of clusters £ and/or the number of Estr used in the 
fit. 

Fitting was performed with the Alloy Theoretic Automated Toolkit (ATAT)^^*^^"^ which 
automates most of the tasks associated with the construction of a CE Hamiltonian. A 
complete description of the algorithms underlying the code can be found in^^. The zero- 
and point-cluster values were -0.421118 eV and 0.006221 eV, respectively. The six pair and 
six 3-body ECI that comprise the complete CE Hamiltonian are plotted in Figs. [2^ and [2)3, 
respectively. ECI for the isostructural TiOx (open symbols, blue online) and HfOx (open 
symbols, red online) systems are also plotted for comparison. As expected, nearest neighbor 
(nn) 0-0 pairs are highly energetic, and therefore strongly avoided; hence nn-pair ECI are 
strongly attractive (ECI >0, for 0-[ ] nn pairs); but beyond nn-pairs, the 0-[ ] pairwise 
ECI are close to zero. The ratio of magnitudes for nn-pair ECIs that are parallel- (Jy) and 
perpendicular (J_l) to cuex, respectively, is J\\/J± ~ 2.5. Note that the 4'th nn-pair ECI is 
the second-nn parallel to chsx, (Jparaiiei) and Jparaiieil J\\ ~ 0.09. 

These results are similar to those presented in Ruban et al.'^ although their effective pair 
interactions and ours are not identically defined. 

Long-period superstructure (LPSS) phases were reported^^ in samples with with bulk 
compositions close to ZrsO (essentially the a'g field in Aral and HirabayashP^, their Fig. 
9). Hence, it is reasonable to speculate that the CE-Hamiltonian might be like that in an 
axial next nearest neighbor Ising model (ANNNI-model),'^in which one expecs Jy and Jl 
to be opposite in sign, and of comparable magnitudes (0.3 ^ — JU J\\ ~ 0.?^; however, 
JUJ\\ ^ 0.09 (Figs. [21). 
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FIG. 2: Effective Cluster Interactions (ECI) for pair and 3-body interactions. Solid black symbols 
indicate the ZrOx-fit which was used in the phase diagram calculation reported here. Open squares 
and down-pointing triangles (blue online) indicate the results of a fit for the TiOx system. Open 
circles and open up-pointing triangles (red online) are from a fit for the HfOx system, (a) The first 
two pair-ECI are for nearest-neighbor 0-[ ] pairs that are parallel- (Jy) and perpendicular (J_l), 
respectively, to Cffex^ and the 4'th nn pair-ECI is the second-nn parallel to Cfjex {J\\)- Pairwise-ECI 
are plotted as functions of inter-site separation, (b) 3-body interactions are plotted as functions 
of the index nj_j_fc which increases, nonlinearly, as the area of triangle i-j-k increases. Positive 
pairwise ECI imply a strong nn-pairwise 0-[ ] attraction, i.e. a strong nn-pairwise 0-0 repulsion. 

III. RESULTS 

A. Ground-States 



The CE was used for a ground-state (GS) analysis that included all configurations of 
[ ] and O in systems of 18 or fewer Zr-atoms (octahedral interstitial sites); a total of 
2^® = 262,144 structures (reduced by symmetry). Five GS were identified in the range, 
< X < 1/2, i.e. at X = 0, 1/6, 1/3, 5/12 and 1/2; sohd circles (black online) on the 



convex hull (solid line) in Fig. [1} The extension of the convex hull towards monoclinic 
zirconia {Zr02) is also plotted in Fig. [I] The CE-results suggest that all four VASP-GS 
in the aZr[ ]i^xOx subsystem are also GS of the Zr-0 binary, and that the maximum 
solubility of O in aZvhcp is X^aa; ~ 1/2 (higher than the experimental value, X ;^ 0.4). 
Note that, the predicted CE-GS at Zr302 is not a GS for the VASP calculations (not a 
VASP-GS); hence the VASP-predicted maximum solubility of O in Zr is Xmax ~ 0.5. 

The larger open squares (red online) in Figure [I] are CE-calculated values for the 
AEf that correspond to the VASP calculations, and the smaller open squares (blue on- 
line) are AEf for the remaining 262,144-96=262048 structures in the GS analysis. The 
open diamond symbol (green online) indicates the calculated formation energy for the P6322 
structure for Zt^O that was originally proposed by Holmberg and Dagerhamn^; this struc- 
ture is also described in Table I. All space group determinations were performed with the 
FINDSYM program.lISMl 



TABLE I: Crystal structure parameters for predicted ground- 
state phases in the aZr[ ]i^xOx system. Cell constants are 
given in A. 



System 


X 


Space Group 


Calculated cell 


Idealized 




atomic 


IT number 


constants 


Atomic 




fraction 


Pearson Symbol 


{A) 


Coordinates 


ZreO 


1/6 


R3 


a ^ \/3ao 


0: 1/6, 1/6, 1/6 






148 


= 5.5333 


Zr: 3/4, 1/12, 5/12 




1/7 


hP7 


c PS 3co = 15.333 


Zr: 11/12, 7/12, 1/4 
Zr: 1/12, 5/12, 3/4 
Zr: 1/4, 11/12, 7/12 
Zr: 5/12, 3/4, 1/12 
Zr: 7/12, 1/4, 11/12 


ZraO 


1/3 


R3c 


a ss y/Soo 


0: 1/6, 1/6, 1/6 






167 


= 5.5671 


0: 2/3, 2/3, 2/3 




1/4 


hP8 


c « 3co = 15.381 


Zr: 3/4, 1/12, 5/12 
Zr: 11/12, 7/12, 1/4 
Zr: 1/12, 5/12, 3/4 
Zr: 1/4, 11/12, 7/12 
Zr: 5/12, 3/4, 1/12 
Zr: 7/12, 1/4, 11/12 


ZrgO 


1/3 


P6322 


a ^ ^/3ao 


0: 1/3, 2/3, 






182 


= 5.5585 


0: 2/3, 1/3, 1/2 




1/4 


hP8 


c« Co = 5.1327 


Zr: 1/3, 0, 
Zr: 0, 1/3, 
Zr: 2/3, 2/3, 
Zr: 2/3, 0, 1/2 
Zr: 0, 2/3, 1/2 
Zr: 1/3, 1/3, 1/2 



Zri205 


5/12 


R3 


a ss y/3ao 


0: 1/12, 1/12, 1/12 






148 


= 5.5568 


0: 1/4, 1/4, 1/4 




5/17 


hP17 


c PS 3co = 30.861 


0: 1/2, 1/2, 1/2 
0: 2/3, 2/3, 2/3 
0: 11/12, 11/12, 11/12 
Zr: 1/8, 11/24, 19/24 
Zr: 1/24, 17/24, 3/8 
Zr: 23/24, 7/24, 5/8 
Zr: 21/24, 13/24, 5/24 
Zr: 19/24, 1/8, 11/24 
Zr: 17/24, 3/8, 1/24 
Zr: 5/8, 23/24, 7/24 
Zr: 13/24, 5/24, 21/24 
Zr: 11/24, 19/24, 1/8 
Zr: 3/8, 1/24, 17/24 
Zr: 7/24, 5/8, 23/24 
Zr: 5/24, 7/8, 13/24 


ZrsO 


1/2 


P31m 


a « a/Soo 


0, 0, 






162 


= 5.5501 


0: 1/3, 2/3, 1/2 




1/3 


hP9 


c« Co = 5.1572 


0: 2/3, 1/3, 1/2 
Zr: 0, 1/3, 3/4 
Zr: 1/3, 1/3, 1/4 
Zr: 1/3, 0, 3/4 
Zr: 2/3, 0, 1/4 
Zr: 2/3, 2/3, 3/4 
Zr: 0, 2/3, 1/4 
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FIG. 3: Idealized crystal structures of the four cluster-expansion-predicted suboxide ground-states: 
(a) ZrgO; (b) ZrsO; (c) Zri205; (d) Zr20. Spheres connected by bond-sticks (yellowish-green online) 
represent Zr. Isolated spheres with bond-sticks (blue online) represent oxygen. Isolated spheres 
(red online) represent vacant octahedral sites. 
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Crystal structures of the VASP- and CE-GS in Zr-ZrO are described in Table I and their 
idealized structures are drawn in Figures [3] a-d: Zr is represented by spheres connected with 
bond-sticks (yellowish-green online); O is represented by isolated spheres with bond-sticks 
(blue online); and [ ] are represented by isolated spheres (red online). 

Various low-T ordered structures have been reported,'^'^^^ with the most recent review 
by Sugizaki et al.^ who used neutron diffraction to study short-range order in ZrOo.s solid 
solutions. Their Figs, la-c presented representations of three ordering schemes that were 
observed within different homogeneity ranges: (a) ZrO^; at X ^ 1/3 (P321); (b) ZrO^, at 
1/3 < X < 0.4 (P6322); (c) ZrO^ near the solubility limit X ^ 0.4 (P31m). Space groups 
for these idealized structures were not reported by Sugizaki et al.™ they were determined in 
this work with the FINDSYM program.'^ Comparing structures (a)-(c) above to the results 
of this work: (a) VASP calculations indicate that this structure is clearly not a GS; (b) is 
the P6322 structure shown as an open diamond (green online) in Fig. [11 its formation energy 
is very close, but higher than, the VASP-GS at X=l/3; (c) is the same P31m structure as 
the VASP-GS at X=l/2. 

1. ZreO, X=l/6, a'i 

The structure of ZrgO is thought to be isomorphic to that of HfeO and TigOB^'^ a~ 
A/3ao; c~co; Z=3 (ao and cq are the cell constants of the disordered PGsmmc alloy) .1^ This 
is also the VASP-GS at X=l/6, Fig. |3](a) and Table 1. 



2. ZrsO, X=l/3, a'^....a'l 

Based on X-ray diffraction studies, Holmberg and DagerhamrP proposed a P6322 struc- 
ture (open diamond, green online, in Fig. \\\ with a?^ y/^do and c~ Cq for a sam- 
ple with X> 0.26. Based on single crystal neutron diffraction studies YamaguchP re- 
ported X-ray, electron and neutron diffraction data on samples in the range ZrO0.i8-ZrO0.30 
(l/5< X < 3/7) and listed atomic coordinates for a "P3cl" structure with a ~ a/Soq, 
c~ 3co. YamaguchP also reported confirmation of the P6322 structure in the composition 
range 0.33 < X <0.4 (1/2 <X<2/3). The FP resuhs presented here suggest that the VASP- 
GS at X=l/3 is the R3c structure depicted in Figure k5l(b). The calculated energy-difference 
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between these two structures is only 0.006 eV, and this difference is probably within DFT 
error, but the precision of these calculations is sufficient to recognize the R3c structure as 
the VASP-GS. 

3. ZruO^, X=5/17 

This structure does not correspond to any reported phase, and because it is predicted to 
disproportionate at T> 20K. It is not expected to be observed experimentally. 

4. ZriO, X=l/2, a'i 

The only Zr20 structure listed in PearsorP^is cubic, and the apparent solubility limit of 
X^ 0.4, rather than X=l/2, which suggests that the VASP calculations may underestimate 
the stability of monoclinic Zr02, and therefore finds the GS tieline between the P31m GS 
at X=l/2 and monoclinic Zr02, rather than between the R3c GS at X=l/3 and monoclinic 
Zr02. Another possibility is that the experimentally measured low-temperature equilibrium 
between Zr-suboxides and monoclinic Zr02 was measured at too low a fugacity of oxygen to 
stabilize the P31m phase at X=l/2. As one expects from the ECI (Fig. ^, there are no 0-0 
nn pairs in the VASP-GS P31m structure, or in any of the four structures with formation 
energies within 0.01 eV (right panel Fig. II]). 

B. The Phase Diagram 

A first principles phase diagram (FPPD) calculation was performed with grand canonical 
Monte Carlo (MC) simulations using the emc2 code which is part of the ATAT package^^"^. 
Input parameters for emc2 were: a simulation box with at least 1568 octahedral sites 
(15x15x6 superceU); 2000 equilibration passes; 2000 Monte Carlo passes. The predicted 
phase diagram is shown in Figure |4j Most phase boundaries were determined by following 
order-parameters of the various ordered phases as functions of X and T; here order param- 
eters are defined such that they are unity in a specified GS-phase, zero in the disordered 
phase, and typically some non-zero value in ordered phases other than their specified GS. 
Dotted boundaries are used to acknowledge uncertainties in phase boundary determinations. 
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FIG. 4: Comparison of experimental and calculated phase diagrams for the system aZr[ ]i^xOx'- 
(a) a combination of the "transformational diagram" (symbols) and the "tentative diagram" (solid 
lines) in Arai and Hirabayashi 1^ (their Figs. 1 and 9, respectively); (b) the diagram calculated from 
this work (dotted phase boundaries are less precisely determined than solid boundaries). Note the 
different results for 0.25 < X < 0.42 and 420iv: < T < 725iC. 



C. The Intermeadiate Temperature a'-Phase 

As observed experimentally in samples with X=0.41,'i^ (up-pointing triangles, green on- 
line, Fig. |4]) a two-step order- disorder process is predicted for 0.25 ^ X ^ 0.5 Figuresjs] The 
data reported in Hirabayashi et al.'^ appear to indicate that both order-disorder transitions 
are second-order (continuous) in character, at least at X=0.41, but the calculations reported 
here suggest that the lower-T transition is strongly first-order (at least at X=l/2) while the 
higher-T transition is continuous. 
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FIG. 5: Calculated order-parameter vs temperature curves for: a) X=0.41; b) X=l/2. Heating 
simulations are indicated by right-pointing triangles (red online) and cooling simulations are rep- 
resented by left-pointing triangles (blue online). As observed experimentally, there is a two-step 
disordering process on heating. 

The simulated intermediate-temperature a'-phase crystal structure was determined 
by symmetry analysis, using the ISODISTORT progranpSEl There are two plausible 
transition paths from the PGs/mmc high-T disordered phase to the P31m GS: 

• (1) PGs/mmc — )■ PGs/mcm — )■ P31m, Ki irrecucible representation, (-1,- 
2,0), (2,1,0), (0,0,1) basis; 

• (2) PGs/mmc — )■ P3ml — )• P31m, P^ irreducible representation, (0,-1,0), (1,1,0), (0,0,1) 
basis. 

Path (1) can be ruled out because it requires a first-order PGs/mmc -^ PGs/mcm transition, 
with unit-cell expansion along both auex axes, which neither experiment nor computation 
supports. 

Path (2) permits a continuous PGs/mmc — )■ P3ml transition, as observed experimentally 
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TABLE II: Atomic positions in PSnil (IT 164) a' crystal 
structue: a ^ ao ^ 3.321; c w cq « 5.14000l; * X=0/Zr. 



Atom 


Wyckoff site 


X 


y 


z 


occupancy 


Zr 


2d 


1/3 


2/3 


..1/4 


1 


Oi 
O2 


la 
lb 











1/2 


xoi<l/2 
2X* - xoi 



and supported computationally. The average a', P3ml structure is described in Table [n] 
and depicted in Fig. [6] where partially occupied 0:[ ]-sites are represented by relatively 
smaller and larger spheres (blue online). The precise occupations of sites Oi and O2 can 
be written as x ^ind 2X — x, respectively; where x < 1/2 is the 0-occupancy of site Oi, 
and X=0/Zr; i.e. at X=0.41 and x = 0.22 then 2X - x = 0.60. With respect to space- 
group determination, the only requirement is that the occupancy of Oi must be different 
from that of O2. The P3ml structure is clearly consistent with the computational results 
shown in Figures [7^ and [TJo. The 0:[ ]-distributions (online 0=red, [ ]=gray) in these 
figures were simulated on reduced (6x6x3) supercells by cooling from lOOOK to 900K. For 
clarity Zr-atoms are omitted to highlight the strong preference for 0:[ ]-ordering along 
CHex] i-G- strong 0-0 nn avoidance along cnex- In the average P3ml structure this leads to 
alternating nn-layers, _L chsx that are relatively 0-rich and 0-poor ([ ]-rich). Visually, this 
statistical difference is obscured in the simulation snapshots (Figures [7^ and [7|d) because 
one has: discrete O and [ ]; 0:[ ]-disorder; and antiphase boundaries. 

IV. DISCUSSION 

A. Comparison of Calculated and Experimental Phase Diagrams 



The main differences between the FPPD presented here and the "tentative phase dia- 
gram" in Aral and HirabayashP (Fig. §1; their Fig. 9) are with respect to their represen- 
tations of broad homogeneity ranges for three ordered phases in the range 0.25 ^ X ^S 0.42 
and 420K ^ T ^ 725K. In this range, Aral and Hirabayashi report three low-T ordered 
phases, a'g, a'^, and a'l; whereas the FPPD has only two; note that the predicted GS at 
X=5/12 disproportionates at T^ 20K. Also, the FPPD-predicted a'-phase field is signifi- 

15 



P3ml a'-structure 




Primative Cell 



FIG. 6: Average PSinl structure of the a'-phase. Small and intermediate sized spheres (blue online) 
represent less- and more oxygen-rich oxygenivacancy-sites (0:[ ]-sites), respectively. Larger spheres 
(red online) represent Zr atoms. More- and less 0-rich 0:[ ]-sites segregate into alternating layers 
perpendicular to chsx', reflecting nearest neighbor 0-0 avoidance. 

cantly larger than the corresponding field in Fig. |4k, and at X=0.41 the a'-aZr^cp transition 
is predicted to occur ^ SOO-ft' higher than experiment suggests, Fig. |5^. Typically, FPPD 
calculations overestimate order-disorder transition temperatures especially when, as here, 
the excess vibrational contribution to the free energy^S! is ignored; so it is not surprising that 
agreement between experiment and theory is not close for the a' ^ aZr^^p order-disorder 
transition. Note however, that the maximum temperatures for stabilities of phases other 
than a' are roughly equal to those shown in Fig. |4^. 



B. Long-Period Superstructures at Xssl/3 

Based on X-ray, neutron, and electron scattering data, Fehlmann et al.'^ and Yamaguchi 
and HirabayashP reported a variety of long-period superstructures (LPSS) in samples with 
bulk compositions X^l/3 (the a'^ field. Fig. |4^) that were subjected to various heat treat- 
ments. The FPPD calculation presented here does not predict LPSS fields, but a similar 
calculation for HfOx appears to predict Devil's Staircases of ordered phases at HfsO and 
Hf20.^° In an ANNI-model like Hamiltonian, one expects, 0.3 ^ —JUJ\\ ^ 0.7, however, the 
12-pair fit which includes Jl yields J\\ and Jl with the same sign and Jl ~ ^||/10. Physically, 
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a)X=0.41,T=900K 



b) X=l/2, T=900K 



FIG. 7: Simulated 0:[ ]-sites (redigray online, respectively) distributions at: (a) X=0.41, T=900K; 
and (b) X=l/2, T=900K. For clarity, Zr-atoms are omitted and a reduced (6x6x3) supercell were 
used. At X=0.41 no 0-0 nn pairs are evident parallel to chsx- At X=l/2, almost all nn pairs 
parallel to chsx are 0-[ ], although two columns (first row, columns 4 and 5) have some 0-0 nn 
pairs), while perpendicular to chcx there are many more 0-0 and [ ]-[ ] nn pairs. 

the fitted values for J\\ and J± are easy to rationalize in terms of 0-0 nn-repulsion, and 
this argues against stable LPSS phases, unless they are stabilized by competition between 
higher-order interactions; e.g. 3'rd and further nn-pair-ECI or multiplet interactions. In 
fact, FPPD calculations for the HfOx system, which has a CE Hamiltonian very similar to 
that for ZrOx, indicate a Devis's Staircase of LPSS phases at HfsO.'^ 

V. CONCLUSIONS 

Ground-State ordered phases are predicted at X=0, 1/6, 1/3, 5/12 and 1/2, but the 
one at X=5/12 is predicted to disproportionate at T:^ 20K, hence it is not expected to be 
observed experimentally. In the range 0.25 ^ X ^ 0.5, in which Aral and HirabayashP^ 
report three phases (a'g, a'^ and a'^) only two are predicted; i.e. the phase fields that Aral 
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and HirabayashP^ draw for a'g and a'^ are predicted to be a single-phase solid solution. 
Figure la clearly indicates that a zeroth order model for octahedral interstitial O: [ ]-ordering 
is one in which first- and second- nn pairwise interactions (Jy and J_l, respectively) strongly 
favor 0-[ ] nn-pairs; i.e. 0-0 nn-pairs are highly unfavorable, and J\\/J± ~ 2.5. Including 
J|l in the ECl fit does not yield an ANNNI-hk^ CE-Hamiltonian; however, recent FPPD 
calculations for the HfOx systeni,'32l (the HfOx-CE is very similar to the ZrOx-CE) predict 
Devis's Staircases of LPSS phases at Hi^O and Hf20. 

The most probable transition path (on cooling) for 0-rich solutions, X ^ 0.4 is FGa/mmc 
— 7- F3ml —7- F31m; hence the average a'-structure has F3ml symmetry with alternating O- 
rich and [ ]-rich layers ± chcx- 
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